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Abstract  We  describe  the  development  of  new  HF  data  assimilation  capabilities  for  our  ionospheric 
inversion  algorithm  called  GPSII  (GPS  Ionospheric  Inversion).  Previously  existing  capabilities  of  this  algorithm 
included  assimilation  of  GPS  total  electron  content  data  as  well  as  assimilation  of  backscatter  ionograms.  In  the 
present  effort  we  concentrated  on  developing  assimilation  tools  for  data  related  to  HF  propagation  channels. 
Measurements  of  propagation  delay,  angle  of  arrival,  and  the  ionosphere-induced  Doppler  from  any  number  of 
known  propagation  links  can  now  be  utilized  by  GPSII.  The  resulting  ionospheric  model  is  consistent  with  all 
assimilated  measurements.  This  means  that  ray  tracing  simulations  of  the  assimilated  propagation  links  are 
guaranteed  to  be  in  agreement  with  measured  data  within  the  errors  of  measurement.  The  key  theoretical 
element  for  assimilating  HF  data  is  the  raypath  response  operator  (RPRO)  which  describes  response  of  raypath 
parameters  to  infinitesimal  variations  of  electron  density  in  the  ionosphere.  We  construct  the  RPRO  out  of  the 
fundamental  solution  of  linearized  ray  tracing  equations  fora  dynamic  magnetoactive  plasma.  We  demonstrate 
performance  and  internal  consistency  of  the  algorithm  using  propagation  delay  data  from  multiple  oblique 
ionograms  (courtesy  of  Defence  Science  and  Technology  Organisation,  Australia)  as  well  as  with  time  series  of 
near-vertical  incidence  sky  wave  data  (courtesy  of  the  Intelligence  Advanced  Research  Projects  Activity  HFGeo 
Program  Government  team).  In  all  cases  GPSII  produces  electron  density  distributions  which  are  smooth  in  space 
and  in  time.  We  simulate  the  assimilated  propagation  links  by  performing  ray  tracing  through  GPSII-produced 
ionosphere  and  observe  that  simulated  data  are  indeed  in  agreement  with  assimilated  measurements. 


1.  Introduction 

HF  sky  wave  propagation  channels  are  known  to  be  greatly  influenced  by  ionospheric  variations.  These 
variations  manifest  themselves  in  HF  measurements  as  fluctuations  of  the  propagation  delay,  angles  of 
arrival,  and  the  ionosphere-induced  shift  of  signal  frequency  (ionospheric  Doppler).  This  sort  of  data  carries 
information  about  spatial  and  temporal  structure  of  the  electron  density  in  the  ionosphere,  albeit  in  a  convoluted, 
implicit  way.  We  have  developed  a  formulation  for  recovering  equivalent  ionospheric  structures  from  measurements 
of  HF  probes. 

This  problem  of  inverting  HF  data  is  approached  as  the  task  of  extending  capabilities  of  GPSII — the 
assimilative  model  of  the  ionosphere  described  in  Fridman  et  al.  [2006,  2012].  GPSII  is  capable  of  ingesting 
data  from  GPS  and  low  Earth  orbiting  satellite  beacons,  in  situ  electron  density,  vertical  incidence  soun¬ 
ders,  and  leading  edges  of  backscatter  ionograms  to  derive  a  three-dimensional  ionosphere  model  that 
is  both  spatially  and  temporally  smooth  but  is  yet  in  agreement  with  all  the  input  data  to  within  the  data 
measurement  error. 

In  the  next  section  we  will  describe  our  formulation  for  incorporating  HF  data  into  ionospheric  inversions 
followed  by  illustrations  of  algorithm  operation.  The  companion  paper  by  Nickisch  et  al.  [2016]  in  this  issue 
offers  quantitative  analysis  of  the  algorithm  performance  for  one  data  collection  campaign. 

2.  Merging  HF  Data  With  GPSII 

Our  approach  to  solving  the  ionospheric  inverse  problem  is  inherited  from  the  GPSII  algorithm  [Fridman 
et  al.,  2006,  2012].  We  represent  the  three-dimensional,  time-varying  distribution  of  electron  density  in 
the  ionosphere  as 


n(r,f)  =  n0(r,  t)Q[u(r,t)} 


(1) 
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where  n0{ r,  t )  is  a  background  model  of  the  ionosphere,  Q(s)  is  a  user-defined  monotonically  increasing  func¬ 
tion  of  one  variable  such  that  Q{—°°)  =  0,  Q(°o)  =  oo,  and  u(r,  t)  is  an  arbitrary  function  which  will  be  determined 
as  a  result  of  the  inversion  procedure.  Results  presented  in  this  paper  employ  IRI-2007  as  the  background 
model.  The  function  Q(s)  is  selected  in  the  following  form: 

{es  ,  s<  0 

1  +s  +  0.5s2-s3/3,  0  <  s<0.5  (2) 

1.25s +  23/24,  s  >  0.5 

This  function  is  continuous  along  with  its  first  and  second  derivatives,  it  decays  exponentially  at  negative  s, 
and  it  behaves  linearly  at  large  positive  s. 

The  numerical  solution  will  be  performed  over  a  four-dimensional  spatial-temporal  grid.  The  vector  of  values 
of  u(r,  t )  in  all  nodes  of  this  grid  will  be  denoted  U( we  are  using  boldface  only  for  vectors  in  three-dimensional 
physical  space).  The  unknown  vector  U  is  related  to  the  vector  of  available  measured  quantities  Y: 

Y  =  M[U }  +  tj.  (3) 

Here  M  is  the  measurement  operator.  It  is  a  known  nonlinear  operator  that  relates  the  ionospheric  model 
(which  is  specified  by  U)  to  theoretical  estimates  of  each  of  the  measured  quantities,  and  vector  ij  represents 
the  noise  of  measurements.  It  is  assumed  that  the  noise  covariance  matrix  of  measurement  errors  is  known. 

The  task  of  the  ionospheric  inverse  problem  is  to  resolve  equation  (3)  with  respect  to  U.  In  order  to  solve  it 
within  the  GPSII  framework,  it  is  necessary  to  provide  L — the  linearized  version  of  the  measurement  operator 
M:  L  =  SM[U]/Su.  Once  procedures  for  calculating  the  vector  M[U]  and  the  linear  operator  L  are  in  place,  the 
unknown  vector  of  ionospheric  modification  U  may  be  found  iteratively  using  Tikhonov's  regularization  tech¬ 
nique  with  the  residual  principle  as  described  in  section  2.2  of  Fridman  et  al.  [2006].  This  inversion  technique 
takes  into  account  the  statistics  of  errors  of  measurements  and  ensures  that  diverse  data  types  are  ingested 
simultaneously  and  with  proper  weights. 

Thus,  in  order  to  incorporate  the  HF  channel  probe  data  into  GPSII,  we  need  to  provide  means  for  calculating 
the  theoretical  channel  probe  data  for  any  ionosphere  represented  by  (1).  This  will  augment  M[U]  with  the 
components  representing  the  HF  data  MHf:[U].  We  also  need  to  augment  L  with  LHf  =  8MHF[U]/Su  which  is 
the  linearization  of  MHf-  Both  tasks  are  solved  with  the  aid  of  numerical  ray  tracing  as  outlined  below. 

2.1.  Extended  Ray  Tracing  Equations 

We  employ  Hamiltonian  formulation  of  the  ray  tracing  (RT)  equations  [Flaselgrove,  1 957;  Jones  and  Stephenson,  1 975]: 

dR  _  8H  dH  dk  _  8H  8H  dro  _  8H  8H 

dr  dk  dco  ’  dr  3R  8co  ’  dr  8t  8co 

where  r  is  the  group  delay,  R(r)  is  the  position  of  the  wave  packet  in  space,  k(r)  is  the  wave  vector,  and  Fi{ R,  k,  ai,  t) 
is  the  RT  Hamiltonian.  Initial  values  of  solution  components  [R,  k,  to]|r  =  0  must  satisfy  the  dispersion  relation 
H(R,  k,  co,  t)  =  0. 

The  Hamiltonian  may  be  in  any  form  discussed  by  Jones  and  Stephenson  [1975].  The  default  form  of  the 
Hamiltonian  in  GPSII  is 


where  //(k,  R)  is  the  refractive  index  of  cold,  collision-free  plasma  as  presented  in  Davies  [1990,  expression  (3.9)]. 

In  order  to  simplify  derivations,  we  are  going  to  use  Cartesian  coordinates  throughout  this  paper. 
Transformation  to  the  Earth-centered  spherical  coordinates  that  are  employed  in  our  practical  codes 
is  straightforward. 

It  is  convenient  to  rewrite  equation  (4)  as  a  single  equation  for  the  seven-component  solution  vector  X  =  [R,  k,  co]T: 


(6) 
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We  explicitly  show  here  that  the  right-hand  side  above  depends  on  spatial  fields  of  electron  density  and  the 
time  derivative  of  electron  density. 

In  order  to  construct  EHF,  we  augment  the  RT  equations  with  their  linearized  version: 

J 

—A  =  BA  (7) 


where  the  seven-element  columns  of  matrix  A(t)  represent  linearized  solutions  of  the  system  and  6  is  a  7  by  7 
matrix  whose  elements  b-,j=dfJ8Xj  are 


hi.: 


-1  T 


8  8H  8H 

- ® —  — 

SR  dkldco 

8  8HI8H 

3R®3R I  dm 

8  8HI8H 

713  “  ^dRetld^’ 


hi  .3,4:6  — 


84:6.4:6  — 


8  8HI8H' 
3k  3k/  3® 

'8  8H I8H 
3k®3R/3/o 
3  8H I8H 


8  8H  I8H 


87,4:6  —  — 


3k  3f  /  3®’ 


^13’7  3®3k/3® 

_  8  8H I8H 
4;6’7  _3®3R/3(u 
8  8H I8H 
7,7  8a>  8t  1 8co 


(8) 


In  our  algorithm  the  extended  system  of  RT  equations  (6)  and  (7)  is  solved  numerically  for  the  7  by  7  matrix  A 
(r)  subject  to  the  initial  condition  A(0)  =  /,  where  /  is  the  identity  matrix.  Thus,  columns  of  A{z)  comprise  the 
fundamental  solution  of  (7). 


2.2.  Raypath  Response  Operator 

For  a  given  receiver-transmitter  configuration  of  a  channel  probe  and  for  a  given  distribution  of  electron  den¬ 
sity  n(r,  f)  one  can  find  rays  that  connect  receiver  and  transmitter  [Coleman,  201 1]  in  accordance  with  the  RT 
equations.  Then  theoretical  values  for  probe's  measurements  at  the  receive  site  (such  as  propagation  delay, 
angles  of  arrival,  and  ionosphere-induced  Doppler  shift)  can  be  expressed  in  terms  of  components  of  the  vec¬ 
tor  J(t)  =  [zL,  kLl  coL],  where  zL  is  the  value  of  z  at  the  landing  point  of  the  traced  ray,  kF  =  k(zL),  and  ®  =  <o(zL). 


The  following  is  the  definition  of  the  raypath  response  operator  (RPRO).  The  RPRO  relates  infinitesimal  varia¬ 
tions  of  n(r,  t)  at  the  time  instant  f  to  variation  of  components  of  the  vector  J\ 


SJ  = 


Rs{r)Sn  +  RD(r)S(-n 


d3r 


(9) 


Plere  Rs  and  R°  are  static  and  dynamic  components  of  the  RPRO. 


The  relationship  between  iHF  and  RPRO  is  straightforward  (as  is  illustrated  at  the  end  of  this  section),  so  we 
will  concentrate  on  deriving  components  of  the  RPRO.  The  plan  is  to  construct  a  Green  function  for  the  linear 
equation  (7).  This  will  allow  us  to  relate  perturbation  of  the  RT  solution  to  infinitesimal  variations  of  the  right- 
hand  side  in  (6).  The  variations  of  the  later  can  be  directly  related  to  perturbations  of  electron  density. 


Assuming  that  the  matrix  A(z)  is  available  as  a  result  of  numerical  ray  tracing,  the  Green  function  (the  impulse 
response  function)  for  (7)  regarded  as  the  initial  value  problem  is 


0,  T<T 

A(r)A_1  (r),  r  >  z 


(10) 


Indeed,  this  function  satisfies  the  impulse  response  equation  for  the  linearized  system 

J 

—  G  =  BG  +  IS(t-t)  (11) 

dr  v  ' 

subject  to  zero  initial  conditions  (G,|I  =  o  =  0).  What  we  actually  need  is  the  Green  function  G(z,z')  (that  is,  the 
solution  to  equation  (11))  for  the  boundary  value  problem  where  end  points  of  perturbed  rays  remain  fixed 
[Afanasyev  et  a!.,  1 982].  In  this  case  the  starting  wave  vector  ks  and  the  group  delay  at  landing  zL  are  affected 
by  perturbations  along  the  ray.  Denoting  the  unknown  response  of  ks  as  K'(r')  (this  is  a  3  x  7  matrix  related  to 
the  yet  unknown  Green  function  as  K(z')  =  G(0,z\4.6j  :7])  and  the  response  of  zL  as  TM  (this  is  a  1  x  7  matrix), 
we  can  express  the  Green  function  of  the  point-to-point  ray  perturbation  problem  as 

G(r,r)  =  G1  (r,  z)  +  G'(r,  0)[1;7>4;6]K(r) [1;3  ,;7]  (12) 
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The  second  term  in  this  expression  introduces  raypath  perturbations  associated  with  perturbations  of  the  starting 
wave  vector  while  position  of  the  starting  point  remains  unperturbed.  The  response  functions  /CM  and  TM  are 
determined  from  the  requirement  that  the  perturbation  in  spatial  position  of  the  landing  point  is  zero: 


(rtl  T) [1.-3, 1:7]  +  G  (r/.,  0)[1;3i4;6]/C(r )  j1;3  1;7j  +  X(r/.)[1:3]T(r  —  0 


(13) 


K(t)  =0 


*(0)[,:3] 

Here  X  =  dX /dr,  found  within  the  unperturbed  RT  problem.  The  second  equation  above  reflects  the  require¬ 
ment  that  the  dispersion  relation  H(R,  k,  co,t)  =  0  remains  fulfilled  by  the  perturbed  starting  wave  vector  (we 
exploited  the  first  equation  in  (4)  to  express  this  requirement  in  a  simpler  form).  Solving  the  matrix  equation 
(13)  with  respect  to  KM  and  TM,  we  find 

G  (t/.,  0)[1:3 4:6]  ^(rz.)[i:3,i] 

0 


K(r') 

T(t') 


/[1: 3, 1:7] 


0 


(14) 


*(<>)[,, 1;3] 

Upon  substituting  KM  into  (12),  we  obtain  the  impulse  response  of  state  vector  X as  a  function  of  the  propagation 
delay  i  from  the  ray  starting  point.  The  response  function  for  the  state  vector  at  the  location  of  the  receiver  is 

Gl{t)  =  G'(  tl,t)  +  G,(ri,0)[1:7i4:6]K(T)  +X(Tt)[1;A1]T(r')[11;7]  (15) 

The  last  term  in  the  above  expression  accounts  for  the  perturbation  of  the  propagation  delay. 

The  components  of  the  RPRO  may  be  expressed  in  terms  of  the  impulse  response  functions  GL  and  T: 


*L 

Rs( r)  =  { 

[1 ,1:7] 

0 

_G  (t)  [4.7,1 :7]  . 

n 

R°(  r)  =  { 

f"(T)  [1,1:7] 

0 

G  (T)  [4:7,1 :7]  _ 

SF(X(z),[n,8n/et})  ^ 
Sn 

dF(X(z),[n,8n/8t}) 

S(8n/8t) 


(16) 

(17) 


In  these  expressions  the  functional  derivatives  are  evaluated  assuming  that  n  and  8n/8t  are  independent 
fields.  The  following  two  properties  of  functional  derivatives  should  be  kept  in  mind  when  evaluating 
SF/Sn  in  (16)  and  (17).  Suppose  that  a  functional  tp[ri\  amounts  to  a  simple  function  (<p[n]=<p(n(t))). 
Then  =  <p(n(r))d(r  —  r) ,  where  <5(r  —  r')  is  the  delta  function  in  three-dimensional  space  and  ip'(n) 
=  8<p(n)/8n.  The  following  property  is  useful  for  handling  terms  in  the  ray  tracing  equations  that 
contain  spatial  derivatives  of  electron  density:  =  ^S(t  —  r). 


It  is  justified  to  set  R^:4]  =  0  because  relative  magnitude  of  the  Doppler  shift  of  operating  frequency  is 
negligibly  small  for  ionospheric  propagation  paths.  The  only  nonnegligible  component  of  R°  is  the 
Doppler  shift  component  which  can  be  rewritten  as  follows: 


Hn(R(r),k(r),ro0) 

3H(R(r),k(r),ro0)/3co0 


R(r))dr 


(18) 


Here  tu0  is  the  frequency  of  the  transmitted  signal  and  Hn  =  8FI/8n.  The  later  expression  is  formulated  with  the 
understanding  that  the  Hamiltonian  depends  on  the  temporal  coordinate  entirely  via  its  dependence  on  n(r,  f); 
that  is,  H{ r,  k,  a>,  t)  =  H(n(r,  f),  k,  ai). 


As  an  example,  consider  transformation  of  RPRO  into  EHf  for  the  case  of  a  single  Doppler  measurement 
(assume  radian  per  second  units)  over  a  single  propagation  path.  Formally, 

Z-hf  =  n0(r,  t)Q[u{r,  f)]  [RsS(t  -  t)  +  FtDd'(i  -  f)]  (1 9) 

where  f  is  the  time  instant  of  the  measurement,  Q’(s)  =  dQ(s)/ds,  and  S'(s)  =  d<S(s)/ds.  The  operator  iHF  needs  to 
be  approximated  over  the  spatial-temporal  grid  of  the  solution.  For  this  purpose  delta  function  operators  are 
replaced  with  appropriate  interpolation  over  the  grid  operators,  and  the  integral  in  (1 8)  is  approximated  by  a 
finite  sum  as  it  is  described  in  Appendix  A. 


3.  GPSII  Operation  With  Channel  Probe  Data 

Data  from  a  network  of  near-vertical  propagation  channels  were  provided  to  us  by  the  Intelligence  Advanced 
Research  Projects  Activity  (IARPA)  HFGeo  Program.  Figure  1  shows  locations  of  the  receiver  and  five 
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Plasma  Frequency  (MHz)  at  220  km;  13  Aug  2013, 14:18  UT 


transmitters  operated  in  Florida.  An 
example  of  range-Doppler  time  ser¬ 
ies  collected  on  one  of  the  links  is 
shown  in  Figure  2.  (Note  that  the 
use  of  the  term  "range"  here  refers 
to  "slant  range"  or  group  path  length, 
that  is,  delay  times  the  speed  of  light.) 
GPSII  inversion  was  performed  using 
only  HF  data  from  all  links  with  a  time 
cadence  of  1  min;  time  step  of  the 
solution  was  set  to  1  min  as  well. 
Correlation-scale  parameters  of  the 
pseudocovariance  matrix  [Fridman 
et  at.,  2006]  were  set  to  the  value  of  1° 
for  both  horizontal  directions,  while 
the  vertical  scale  varied  from  about 
25  km  at  the  bottom  boundary  of 
the  domain  (altitude  of  80  km)  to 
about  200  km  at  the  top  boundary 
(1000  km).  Solution  grid  size  was 
approximately  0.1°  in  both  horizontal 
dimensions,  while  the  altitude  step 
size  varied  from  about  0.3  km  to  20  km.  The  contour  plot  in  Figure  1  shows  a  snapshot  of  the  horizontal  cross 
section  of  the  solution  at  altitude  220  km  (approximate  altitude  of  reflection  of  the  rays).  One  can  see  that  the 
solution  remains  quite  smooth  while  the  ionosphere  is  modified  to  match  the  data.  The  GPSII  fit  to  the  data  as 
reported  by  the  internal  ray  tracing  is  also  shown  in  Figure  2.  The  GPSII-derived  prediction  follows  data  quite 
closely  during  most  of  the  analyzed  time  interval.  Statistical  analysis  of  GPSII  performance  for  these  kinds  of 
data  is  provided  in  the  companion  paper  by  Nickisch  et  at.  [2016]. 


86°W 


84°w  82UW  80“W 
Longitude  (Deg) 


78°W 


Figure  1 .  HF  transmitters  (circles)  and  receiver  (crosses)  providing  range-Doppler 
data  over  five  propagation  channels.  The  contour  plot  presents  the  horizontal 
cross  section  (220  km  altitude)  of  one  frame  (1418  UT,  13  August  2013)  of  the 
GPSII  solution  driven  by  the  HF  data  from  the  five  transmitters. 


A  valuable  data  set  from  a  network  of  oblique  ionogram  (Ol)  links  and  vertical  sounders  was  provided  to  us  by 
Defence  Science  and  Technology  Organisation  (DSTO)  (Australia).  Figure  3  illustrates  the  configuration  of  the 
Ol  links,  vertical  sounders,  and  GPS  receivers  that  were  employed  in  one  of  our  tests.  These  selected  links  are  a 
small  subset  of  dozens  Ol  paths  monitored  by  DSTO.  Ol  data  were  provided  to  GPSII  as  a  number  of  delay- 
frequency  samples  for  the  extraordinary  ( X)  mode  for  each  participating  Ol  link.  An  example  of  data  samples 
collected  from  one  of  the  Ols  (ASP-SCFI  link)  is  indicated  by  boxes  in  Figure  4  (left).  The  GPSII  solution  was 
produced  with  a  cadence  of  15  min.  GPS  total  electron  content  (TEC)  data  were  sampled  with  15  min  time 
interval.  Ol  and  VI  data  were  occasionally  decimated  to  ensure  that  any  given  Ol  or  VI  contributes  only  once 
per  solution  time  step.  Correlation-scale  parameters  of  the  pseudocovariance  matrix  were  set  to  5°  for  both 
horizontal  dimensions.  The  altitude-dependent  vertical  correlation  scale  varied  from  about  25  km  at  the 


Malabar,  13-Aug-2013  Malabar,  13-Aug-2013 


Figure  2.  Example  of  (left)  range  and  (right)  Doppler  data  compared  to  GPSII  fit  (circles). 


FRIDMAN  ET  AL. 


ASSIMILATION  OF  HF  MEASUREMENTS 


5 


©AGU  Radio  Science 


10.1 002/201 5RS005890 


Synthetic  Ols  were  generated  by  ray  tra¬ 
cing  through  GPSII  ionosphere  using  an 
external  numerical  ray  tracing  code.  A 
good  match  to  the  observed  01  trace  is 
evident  in  Figure  4  when  Ol  data  are 
assimilated  (yellow  points).  It  is  notable 
that  the  same  GPSII  model  provides 
good  agreement  for  the  nonassimilated 
01  link  (WYN-LYN).  Figure  4  also  illus¬ 
trates  that  the  ionospheric  model 
obtained  without  any  01  data  assimi¬ 
lated  (red  points)  fails  to  provide  satis¬ 
factory  agreement  with  the  01  data. 
(Note  that  GPS  L1/L2  beacon  data  and 
vertical  ionograms  were  assimilated  in  all 
cases.  We  also  note  that  the  ionosphere 
was  exceptionally  dynamic  at  this  time  due  to  large-scale  traveling  ionospheric  disturbance  activity.)  Figure  5 
shows  comparison  of  zonal  distributions  of  plasma  obtained  with  and  without  01  data  for  the  time  of  observations 
in  Figure  4.  Figure  5  (bottom)  represents  inversion  without  01  data.  One  can  note  prominent  zonal  modulation  of  the 


Figure  3.  Assimilated  01  links  (red  segments),  assimilated  vertical  profiles 
(triangles),  and  GPS  TEC  receivers  along  with  the  test  01  link  (WYN-LYN). 


bottom  boundary  of  the  domain  (alti¬ 
tude  of  80  km)  to  about  1 000  km  at  the 
top  boundary  (20000  km).  The  solution 
grid  size  was  approximately  0.5°  in  both 
horizontal  dimensions,  while  the  alti¬ 
tude  step  size  varied  from  about 
0.3  km  to  200  km. 


Figure  4.  Ols  collected  on  16  April  2008,  0202  UT,  compared  to  X-mode  synthetic  traces  (yellow  dots),  (left)  One  of  the  assimilated  data  links.  Frequency-delay  data 
supplied  to  GPSII  are  marked  by  green  boxes,  (right)  Data  from  a  validation  (nonassimilated)  link.  Red  dots  show  synthetic  Ols  for  a  different  GPSII  solution  that  was 
obtained  without  assimilating  Ol  data  (only  GPS  and  VI  data  were  assimilated  in  this  case). 
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Plasma  freq.  (MHz)  at  Latitude  -18.00;  16  Apr  2008.  02:02  UT 
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Figure  S.  Vertical  cross  section  of  GPSII  solutions  at  a  constant  latitude  of  1 8°  south,  (top)  Solution  driven  by  Ol,  VI,  and  GPS 
TEC  data  compared  to  the  (bottom)  solution  driven  by  VI  and  GPS  TEC. 


solution  that  may  be  attributed  to  the  impact  of  data  from  nearby  VI  sounders  (at  longitudes  1 1 8.6, 1 23.8, 1 30.8,  and 
1 44.9°)  as  well  as  impact  from  the  line  of  sight  of  one  of  the  GPS  receivers.  The  scale  of  this  modulation  is  consistent 
with  the  adopted  solution  correlation  scale  of  5°.  This  modulation  effect  becomes  less  prominent  with  the  addition 
of  Ol  data  (Figure  5,  top).  It  appears  that  Ol  data  that  probe  geographically  extended  volumes  of  the  ionosphere  help 
to  smoothly  connect  localized  effects  of  sparse  VI  sounders  and  GPS  receivers. 

4.  Conclusions 

We  have  developed  theoretical  framework  for  incorporating  HF  channel  probe  data  (propagation  delay,  angles  of 
arrival,  and  Doppler  shift)  into  ionospheric  inversion  algorithms.  We  have  extended  capabilities  of  GPSII  (our  iono¬ 
spheric  inversion  algorithm)  to  incorporate  data  from  FIF  channel  probes  and  oblique  ionograms.  Performance  of 
the  new  algorithm  was  demonstrated  using  range-Doppler  time  series  produced  by  a  network  of  short-range  FIF 
channel  probes.  In  another  demonstration  we  operated  the  algorithm  with  data  from  a  network  of  Ol  sounders. 
We  performed  ray  tracing  through  the  GPSII-produced  ionosphere  using  an  independent  numerical  ray  tracing 
code  and  observed  that  simulated  data  are  indeed  in  agreement  with  assimilated  measurements. 

Appendix  A:  Numerical  Approximation  of  Delta  Function  Operators  and  Matrix 
Representation  of  LHf 

The  one-dimensional  delta  function  operator  is  defined  by  requiring  that  the  following  relationship  is  valid 
for  an  arbitrary  continuous  function  f(s ): 

\-J{s~p)f{s)ds  =  f{p)  (A1) 

Assuming  that  the  function  f(s)  is  specified  by  a  vector  of  its  values  fj  =  f(s,)  in  nodes  s„  the  operator  (A1 )  may 
be  approximated  by  the  sum 


L*  -  P)f(s)dS*V  ( Sj,rn  Wl  P  +  <W  P  Sm  )  f, 

j\  sm+ 1  —  sm  Sm+1  —SmJ 


(A2) 
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Where  <5 ,_m  is  the  Kroneker  delta,  m  =  m(p )  specifies  the  interval  that  contains  point  p,  and  m :  sm  —P<  sm  + 1  ■  It 
is  easy  to  see  that  (A2)  provides  the  linear  interpolation  approximation  to  f[p).  A  sampled  approximation  to 
the  operator  S’{s  —  p)  is  obtained  by  differentiating  (A2)  with  respect  to  p: 


\_J(s-p)f(s)ds*J2 


dj.m  di  m  I  1 
$m+ 1  $m 


fl 


(A3) 


Multidimensional  delta  function  operators  and  their  first  derivatives  are  obtained  by  applying  transforma¬ 
tions  (A2)  and  (A3)  to  each  of  the  spatial-temporal  dimensions.  In  order  to  illustrate  this  procedure,  we  will 
proceed  with  the  example  of  a  single  Doppler  measurement  that  resulted  in  the  expression  (18)  for  RPRO. 
Substitute  (18)  into  (19)  (keep  in  mind  that  Rs  =  0  in  this  case),  change  order  of  integration  operations,  approx¬ 
imate  integration  over  z  (along  the  raypath)  by  a  sum  over  variable  /  using  the  trapezoid  formula,  and  exploit 
approximations  (A2)  and  (A3).  We  obtain 


(tN+ 1  ,  ,  , 

Jt#l  df J]Jdr/.Hf (r,  f )<5u(r,  f )« 


J2DlJlJlJlJlno(Xii+h^ii+g’Zki+e’  fw+p)Q'  \u(xii+h,y]l+g,Zkl+e,  fN+P)|  M)P+1 

I  p= 0  h= 0  g= 0  e=0 


K+1  -h  -  Rx{l)  |  KJi+1-9  |  kfe,+1-e  - 

xi/+1  —  xh  Yji+ 1  ~  y/,  zh+1  —  zk. 


<5u('/  +  h,ji  +g,ki  +  e,N  +  p) 


(A4) 


Here  [fw,fN+1]  is  the  solution  time  interval  that  contains  the  time  instant  f  of  the  measurement,  [Rx(l),Ry(l),Rz 
(l)]T  =  R(r/),  and  (/,,;,,  kj)  specifies  position  of  the  voxel  of  three-dimensional  grid  of  the  spatial  variable  r=[x,y,z]r 
that  contains  the  point  R (zj),  so  that 


Xi<Rx(l)  <  x/(+, ,  yj <Ry(t)  <  yJ(+1 ,  </?z(/)  <  zk)+1 

D=-GL(z)  Hn(R(T/)»k(r<):w o)  ^<+1  -  r/-1  1 

r/  [7’7l0H(R(r/),k(r,),coo)/3cuo  2  f/v+i  -  f/v 


(A5) 

(A6) 


Expression  (A4)  provides  the  desired  matrix  representation  of  the  operator  /.Hf-  This  matrix  specifies  linear 
response  of  measured  Doppler  to  variations  of  the  function  u  in  nodes  of  the  four-dimensional  grid  of 
the  solution. 
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